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Abstract 



Using a simple transfer matrix approach we have derived very long series expan- 
sions for the perimeter generating function of three-choice polygons. We find that 
all the terms in the generating function can be reproduced from a linear Fuchsian 
CN ' differential equation of order 8. We perform an analysis of the properties of the 

^ ■ differential equation. 

m 

^ '. 1 Introduction 
in 

A well-known long standing problem in combinatorics and statistical mechanics is to 
find the generating function for self-avoiding polygons (or walks) on a two-dimensional 
lattice, enumerated by perimeter. Recently, we have gained a greater understanding of the 
difficulty of this problem, as Rechnitzer 13J has proved that the (anisotropic) generating 
function for square lattice self-avoiding polygons is not differentiably finite |14j . as had 
; been previously conjectured, on numerical grounds but not proved. That is to say, it 

cannot be expressed as the solution of an ordinary differential equation with polynomial 
coefficients. There are many simplifications of this problem that are solvable [Tj, but all 
the simpler models impose an effective directedness or other constraint that reduces the 
problem, in essence, to a one-dimensional problem. 

One model, that of so-called three-choice polygons, has remained unsolved despite the 
knowledge that its solution must be D-finite. In this paper we report on recent numerical 
work resulting in an exact differential equation apparently satisfied by the perimeter gen- 
erating function of three-choice polygons. While our results do not constitute a rigorous 
mathematical proof the numerical evidence is compelling. 

Three-choice self-avoiding walks on the square lattice, Z 2 , were introduced by Manna 
and can be defined as follows: Starting from the origin one can step in any direction; 
after a step upward or downward one can head in any direction (except backward); after a 
step to the left one can only step forward or head downward, and similarly after a step to 
the right one can continue forward or turn upward. Alternatively put, one cannot make a 
right-hand turn after a horizontal step. Whittington jT^j showed that the growth constant 
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Staircase 



Imperfect staircase 



Figure 1: Examples of the two types of three-choice polygons. In the right panel we 
indicate the origin (O) and the direction of the first step (note that rotation by 180 
degrees also leads to a valid three-choice polygon). 

for three-choice walks is exactly 2, so that if w n denotes the number of such walks of n 
steps on an infinite lattice, equivalent up to a translation, then w n ~ 2 n+o( - n \ It is perhaps 
surprising that the best known result for the sub-dominant term is 2°^ but attempts 
to improve on this have not been successful. Even numerically, there is no firmly based 
conjecture for the sub-dominant term, unlike for ordinary self-avoiding walks, for which 
the sub-dominant term is widely believed to be 0(n 9 ). 

As usual one can define a polygon version of the walk model by requiring the walk to 
return to the origin. So a three-choice polygon [Hj is simply a three-choice self-avoiding 
walk which returns to the origin, but has no other self-intersections. There are two 
distinct classes of three-choice polygons. The three-choice rule either leads to staircase 
polygons or imperfect staircase polygons [3] as illustrated in figure ^ In the case of 
staircase polygons any vertex on the perimeter can act as the origin of the three-choice 
walk (which then proceeds counter-clockwise), while for imperfect staircase polygons there 
is only one possible origin but the polygon could be rotated by 180 degrees. If we denote 
by t n the number of three-choice polygons with perimeter 2n then, t n = 2nc n + 2p n , 
where c n is the number of staircase polygons with perimeter 2n, and p n is the number of 
imperfect staircase polygons of perimeter 2n. Note that t n , p n and c n all grow like 4 n and 
in particular we recall the well-known result that c n+ \ = C n = ^r[( 2 ™) are given by the 
Catalan numbers C n . 

In this paper we report on recent work which has led to an exact Fuchsian [§| linear 
differential equation of order 8 apparently satisfied by the perimeter generating function, 
T(x) = ^2 n>0 t n x n , for three-choice polygons (that is T(x) is conjectured to be one of the 
solutions of the ODE, expanded around the origin). The first few terms in the generating 
function are 

T(x) = Ax 2 + 12a; 3 + 42x 4 + 152a; 5 + 562a; 6 + • • ■ . 

The generating function for the coefficients p n is no simpler. 

If we distinguish between steps in the x and y direction, and let £ m>n denote the 
number of three-choice polygons with 2m horizontal steps and 2n vertical steps, then the 
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anisotropic generating function for T can be written 



T(x,y) = Y,^nX m y n = Y,H n (x)y n , 

m,n n 

where H n (x) = is the (rational [T^]) generating function for three-choice polygons 

with 2n vertical steps. In earlier, unpublished, numerical work, we found that, for imper- 
fect staircase polygons, the denominators were: 

S n ( x ) = (1 - xf n -\l + x) {2n ~ 7)+ n even, 

and 

S n ( X ) = (1 - X) 2n -\l + x) (2 "- 8)+ 71 odd. 

This was subsequently proved by Bousquet-Melou [2|. Further, Bousquet-Melou showed 
that the numerators satisfied: 

, > -12(4m)! 

Rn(-i) = — — '— n = 2m + 4, 

nK ' m!(m + l)!(m + 2)!(m + 3)! ' 



and 



, > -96 4m+l! 

RJ-l) = — — ) n = 2m + 5. 

v ; m! m + 1 ! m + 2 ! m + 4 ! 



Unfortunately, we still do not have enough information to identify the numerators, though 
we observe that they are of degree 3n — 7 for n > 4 and n even, and of degree 3n — 8 for 
n > 5 and n odd. 

It is also possible to express the generating function T(x) as a five-fold sum, with 
one constraint, [2. of 4 x 4 Gessel-Viennot determinants [3]. This is clear from figure El 
where the enumeration of the lattice paths between the dotted lines is just the classical 
problem of 4 non-intersecting walkers, and these must be joined to two non-intersecting 
walkers to the left, and to two non-intersecting walkers to the right. Then one must sum 
over different possible geometries. However the fact that the generating function is so 
expressible implies that it is differentiably finite [TT] . 

In the following sections we report on our work leading to an ODE for the perimeter 
generating function of three-choice polygons. We started by generating the counts for 
three-choice polygons up to half-perimeter 260. Using numerical experimentation we 
then found what we believe to be the underlying ODE. This calculation required the use 
of the first 206 coefficients with the resulting ODE then correctly predicting the next 54 
coefficients. While the possibility that this ODE is not the correct one is extraordinarily 
small, our result does of course not constitute a proof. Unfortunately we cannot usefully 
bound the size of the underlying ODE, otherwise we could use the knowledge of D- 
finiteness to provide a proof of our results. That is to say, any bounds that follow from 
closure theorems [TT] are too large to be useful. 



2 Computer enumeration 

The algorithm we use to count the number of imperfect polygons is a slightly modified 
version of the algorithm of Conway et al. 3J. Before proceeding to the full problem it is 
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Figure 2: Showing the decomposition of an imperfect staircase polygon into a sequence 
of 2-4-2 non-intersecting walkers, each expressible as a Gessel-Viennot determinant 

useful to briefly outline the transfer matrix algorithm for enumerating staircase polygons. 
Recall that a staircase polygon consists of two directed walks starting at the origin, moving 
only to the right and up, and terminating once the walks join at a vertex. If we look at a 
diagonal line x + y = k + 1/2 then for any integer k this line will intersect a polygon at 
(miss the polygon) or 2 edges (intersect the polygon), see figure HJ We start with k = 
such that the line intersects the first two edges of the staircase polygon. We then move the 
line upward (increase k by 1) and as we do this we add an edge to each walk. There are 
only four new configurations corresponding to the four possible steps. We need only keep 
track of the gap between the two walks, where the gap is the minimal number of iterations 
required in order to join the two walks. As we move the line the gap is either increased 
by a unit (the upper walk moves up and the lower walk moves right), decreased by a unit 
(the upper walk moves right and the lower walk moves up) or remains constant in two 
possible ways (both walks move up or right). Let C(i, k) be the number of configurations 
with a gap of i at step k. We then have the following very simple algorithm: 

Set C(l, 0) = x (where x is a variable conjugate to the half-perimeter of the polygon). 
Run through all possible gaps % = 1, . . . , k + 1 and do the following updates: C{i + 
1, k + 1) := xC(i, k), C(i - 1, k + 1) := xC(i, k) and C(i, k + 1) := 2xC(i, k). Here a := b is 
short-hand for assign to a the value a + b. 

Formally we can view the transformation from the set of states C(i, k) to C(j, k + 1) 
as a matrix multiplication (hence our use of the nomenclature transfer matrix algorithm) 
with k counting the number of iterations of the transfer matrix algorithm. However, as 
can be readily seen from the algorithm the transfer matrix is extremely sparse and there 
is no reason to list it explicitly (it is given implicitly by the updating rules). 

The term C(0,k) is the number of staircase polygons of half-perimeter k+1. Note that 
the use of the variable x is somewhat superfluous in the case of staircase polygons since 
the generating function at iteration k is just x k+1 C(i, k), where C(i, k) is the number of 
configuration with gap i after k iterations. But it is included here for reasons of generality 
and in the case of imperfect staircase polygons the generating function will be a (non- 
trivial) polynomial in x. Naturally we need not actually keep all the entries C(i, k) since 
only the current and subsequent values are needed for the calculation so we can replace 
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C(i, k) with C(i, k mod 2). We just have to initially set to zero all entries in the next step 
and keep a running total c(k) of the number of staircase polygons. 

Imperfect staircase polygons start out as ordinary staircase polygons (see figure |2J) 
Then at some vertex two additional directed walks (sharing the same starting point) are 
inserted between the two original walks (at the first dashed line marked with a '4' in 
figure [21). The diagonal line will thus intersect these polygon configurations at four edges. 
Imperfect staircase polygons are created by connecting the first two walks and the last two 
walks (as illustrated at the last dashed line marked with a '4' in figure|21). With four walks 
we need to retain three pieces of information, namely, the three gaps /, m, and n between 
consecutive walks. Each existing configuration can produce 16 new configurations as each 
walk is extended by a step either up or to the right. The resulting updating is easily 
worked out jH]- Let G(l,m,n) be the generating function (a polynomial in the variable 
x) for partially completed polygons at a given diagonal k. As we proceed to the next 
diagonal HI we add x 2 G(l,m,n) (the factor x 2 arise because we extend all walks by a 
step) to G(l, m, n) (twice), G(l + 1, m, n), G(l, m + 1, n), G(l, m, n + 1), G{1 — 1, m, n), 
G(l,m — 1, n), G(l,m,n — 1), G{1 + 1, m — 1, n), G(l + l,m,n — 1), G{1 — 1, m + 1, n), 
G(l,m + l,n - 1), G(l - l,m,n + 1), G(l,m - l,n + 1), G{1 - 1, m + 1, n - 1) and 
G{1 + l,77i — l,n + l). Any update resulting in G(l,0,n) has to be rejected because it 
corresponds to a configuration in which we have joined the two middle walks in and this 
can never lead to an imperfect staircase polygon. Obviously once any two walks have 
been connected the remaining walks follow the usual staircase polygon updating rules. 

The configurations with two walks already connected can also be encoded by the G 
functions. We simply let G(l, 0, 0) be the generating function for partial polygons with 
two walks already connected (note that if the boundary line intersects four edges m > 0). 
So in the updating of imperfect staircase polygons we can set G(0, m, n) (we connect the 
two lower-most walks) to G(n, 0, 0) (this case is illustrated at the last dashed line marked 
with a '4' in figure |2I). Likewise we can set G(l,m, 0) (we connect the two uppermost 
walks) to G(l,0, 0). The condition for the formation of a valid polygon is I = n = 
(note that we can't demand m = as well, since we could connect both the lower- and 
uppermost walks simultaneously). 

The 'creation' of a configuration with three gaps, alternatively, one in which a diagonal 
line intersects four edges of an imperfect staircase polygon, which we refer to as a G-type 
configuration, is also very simple (see the first dashed line marked with a '4' in figure |2J). 
We start with a staircase type configuration C(i) and from this we can create four G- 
type configurations by assigning the value x 2 C{i) to G(j,l,i — j), G(j — 1,1, i — j), 
G(j, 1, i—j — 1) and G(j — 1, 1, i— j — 1), where 1 < j '• < i — 1 (the factor x 2 arises because 
we extend the outer walks by a step and insert two new walks each containing a single 
step). 

The algorithm outlined above is already very efficient, but it can be further enhanced 
by the following simple observation. If we wish to calculate the number of polygons up 
to a given maximal half-perimeter length N, we need not consider all possible gaps since 
some configurations can only lead to polygons of a size exceeding N. First of all since gaps 
only increase or decrease by one at each iteration we need never consider configurations 
with gaps exceeding N/2. Furthermore, any G-type configuration with m > must have 
half-perimeter at least k + m. Here we get the contribution k from the outermost walks 
(k is the number of forward steps or iterations taken) and the contribution m from the 
innermost walks (a gap m requires at least m steps). In order to produce an imperfect 
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staircase polygon we have to add at least I + n additional steps (we have to join both the 
two upper-most and two lower- most walks) , so if M = k + l + m + n>N we can discard 
this configuration. Not only can we thus discard some configurations when M > N 
but we can also further decrease the memory use since rather than storing N terms per 
configuration we only need to store N — M terms. 

We calculated the number of imperfect staircase polygons up to perimeter 520. The 
integer coefficients become very large so the calculation was performed using modular 
arithmetic [TU] . This involves performing the calculation modulo various prime numbers 
Pi and then reconstructing the full integer coefficients at the end. We used primes of the 
form pi = 2 30 — r,, where are small positive integers, less than 1000, chosen so that pi 
is prime, and pi 7^ pj unless i = j. 18 primes were needed to represent the coefficients 
correctly. The calculation for each prime used about 250Mb of memory and about 18 
minutes of CPU time on a 2.8 GHz Xeon processor. Naturally we could have carried the 
calculation much further but as we shall demonstrate in the next section this more than 
sufficed to identify an exact differential equation satisfied by T(x). 

3 The Fuchsian differential equation 

In recent papers Zenine et al. [T71 HBJ EH] obtained the linear differential equations whose 
solutions give the 3- and 4-particle contributions an d X^ 4 "* to the Ising model suscep- 
tibility. In this paper we use their method to find a linear differential equation which 
has as a solution the generating function T(x) for three-choice polygons. We briefly 
outline the method here. Let us assume we have a function F(x) with a singularity at 
x — x c = Starting from a (long) series expansion for the function F(x) we look for 
a linear differential equation of order m of the form 



such that F(x) is a solution to this homogeneous linear differential equation, where the 
Pk{x) are polynomials. In order to make it as simple as possible we start by searching 
for a Fuchsian jU] equation. Such equations have only regular singular points. There 
are several reasons for searching for a Fuchsian equation, rather than a more general 
D-finite equation. Computationally the Fuchsian assumption simplifies the search for a 
solution. One may also argue, less precisely, that for "sensible" combinatorial models 
one would expect Fuchsian equations, as irregular singular points are characterized by 
explosive, super-exponential behaviour. Such behaviour is not normally characteristic of 
combinatorial problems. (The point at infinity may be an exception to this somewhat 
imprecise observation). One may also ask the question whether most of the problems in 
combinatorics with D-finite solutions have Fuchsian solutions? While we have not made 
an exhaustive study, we know of no counter-example to this suggestion. 

From the general theory of Fuchsian jU] equations it follows that the degree of Pk{x) 
is at most n — m + k where n is the degree of P m (x). To simplify matters (reduce the 
order of the unknown polynomials) it is often advantageous to explicitly assume that 
the origin and x = x c are regular singular points and to set Pk(x) = Qk(x)S(x) k , where 
S(x) = xR(x) and R(x) is a polynomial of minimal degree having x c as a root (in our 
case we have R(x) = 1 — Ax). S(x) could be generalised to include more regular singular 




(1) 
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points if some were known from other methods of analysis, but we have not found this 
to be advantageous. Thus when searching for a solution of Fuchsian type there are only 
two parameters: namely the order m of the differential equation and the degree q m of 
the polynomial Q m (x). Let p be the degree of S(x) (2 in our case), then for given 
m and q m there are L = (m + l)(q m + 1) + pm(rn + l)/2 — 1 unknown coefficients, 
where we have assumed without loss of generality that the leading order coefficient in 
P m {x) = Qm(x)S(x) m is 1. We can then search systematically for solutions by varying m 
and q m . In this way we first found a solution with m = 10 and q m = 12, which required the 
determination of L = 206 unknown coefficients. We have 260 terms in the half-perimeter 
series and thus have more than 50 additional terms with which to check the correctness 
of our solution. Having found this conjectured solution we then turned the ODE into a 
recurrence relation and used this to generate more series terms in order to search for a 
lower order Fuchsian equation. The lowest order equation we found was eighth order and 
with q m = 30, which requires the determination of L = 321 unknown coefficients. Thus 
from our original 260 term series we could not have found this 8 th order solution since 
we did not have enough terms to determine all the unknown coefficients in the ODE. 
This raises the question as to whether perhaps there is an ODE of lower order than 8 
that generates the coefficients? The short answer to this is no. Further study [7j of our 
differential operator revealed that it can be factorised. In fact we found a factorization 
into three first-order linear operators, a second order and a third order. The generating 
function is a solution of the 8 th order operator, not of any of the smaller factors. 

So the (half)-perimeter generating function T(x) for three-choice polygons is conjec- 
tured to be a solution of the linear differential equation of order 8 

X>(*)|^(z) = (2) 

k=0 

with 



P 8 (x) = x 3 (l - 4x) 4 (l + 4x)(l + Ax 2 ){l + x + 7x 2 )Q 8 (x), 

p 7 (x) = x 2 (l - 4x) 3 Q 7 (x), P 6 (x) = 2x(l - 4x) 2 Q 6 (x), 

P B (x) = 6(1 - 4x)Qb(x), P^x) = 24Q 4 (x), (3) 

P 3 (x) = 2AQ 3 {x), P 2 {x) = 144x(l - 2x)Q 2 {x), 

P x {x) = 144(1 - Ax)Qi(x), P (x) = 576Q {x), 

where Qs(x), Qr(x), . . ., Qo(x), are polynomials of degree 25, 31, 32, 33, 33, 32, 29, 29, 
and 29, respectively. The polynomials are listed in Appendix A (note that the polynomials 
do not factorise). 

The singular points of the differential equation are given by the roots of Pg(x). One 
can easily check that all the singularities (including x = oo) are regular singular points 
so equation (J2J is indeed of the Fuchsian type. It is thus possible, using the method 
of Frobenius, to obtain from the indicial equation the critical exponents at the singular 
points. These are listed in Tabled 

For equations of the Fuchsian type the critical exponents satisfy a simple Fuchsian 
summation relation, which we now take the opportunity to confirm in our case. Let 
Xx, X2, ■ ■ -x n , x n+ i = oo be the regular singular points of a Fuchsian type equation of 
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Table 1: Critical exponents for the regular singular points of the Fuchsian differential 
equation satisfied by T(x). 



Singularity 



Exponents 



x = 



-1, 0, 0, 0, 1, 2, 3, 4 

-1/2, -1/2, 0, 1/2, 1, 3/2, 2, 3 

0, 1, 2, 3, 4, 5, 6, 13/2 

0, 1, 2, 3, 4, 5, 6, 13/2 



x = 1/4 
x = -1/4 
x = ±i/2 
l + x + 7x 2 



0, 1, 2, 2, 3, 4, 5, 6 

-2, -3/2, -1, -1, -1/2, 1/2, 3/2, 5/2 
0, 1, 2, 3, 4, 5, 6, 8 



X = oo 



Q 8 (z) = 



order m and a^i, . . . aj, m (j = 1, . . . n + 1) the m exponents determined from the roots 
of the indicial equation for each regular singular point, Xj, then the following Fuchsian 
relation holds: 



In this case the number of regular singular points is m + 1 = 33, namely the 25 roots of 
Q 8 (x), the two roots of 1 + x + 7x 2 , x — ± i/2, x — ±1/4, x — and x = oo. It is easy to 
verify that the Fuchsian relation is satisfied with m = 8, n = 32, and all the exponents 
aj t k summing to 868, which is a useful check on our results. 

We shall now consider the local solutions of the differential equation around each 
singularity. Recall that in general it is known [3J |U] that if the indicial equation yields 
k critical exponents which differ by an integer, then the local solutions may contain 
logarithmic terms up to log^" 1 . However, for the Fuchsian equation (0) only multiple 
roots of the indicial equation give rise to logarithmic terms in the local solution around 
a given singularity, so that a root of multiplicity k gives rise to logarithmic term up to 



In particular this means that near any of the 25 roots of Qs(x) the local solutions 
have no logarithmic terms and the solutions are thus analytic since all the exponents are 
positive integers. The roots of Qg(x) are thus apparent singularities jHH] of the Fuchsian 
equation (J2J). There are methods for distinguishing real and apparent singularities (see, 
e.g, |3] §45) and in principle one should check that the roots of Qs(x) satisfy the conditions 
for being apparent singularities. However, this theoretical method is quite cumbersome. 
An easier numerical way to see that the roots of Qg(x) must be apparent singularities is 
as follows. We already found a 10th order Fuchsian equation for which the polynomial 
Pio(x) was of a form similar to Ps(x) as listed in equation Q, but with the degree of 
Qio(x) being only 7. That is all the singularities as tabulated in Tabled also appear in 
this higher order equation with the exception of the 25 roots of Qs(x) (at most 7 of these 
could appear in the order 10 Fuchsian equation). In fact we can find a solution of order 14 
of the same form as above but with Qu(x) being just a constant. So at this order none of 
the roots of Qs(x) appear. Clearly any real singularity of the system cannot be made to 
vanish and we conclude that the 25 roots of Qs(x) must indeed be apparent singularities. 

Assuming that only repeated roots give rise to logarithmic terms, and thus that a 




(4) 



j=i k=i 
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sequence of positive integers give rise to analytic terms, then near the physical critical 
point x = x c = 1/4 we expect the singular behaviour 

T(x) ~ A(x)(l - 4x)- 1/2 + B{x){\ - 4x)~ 1/2 log(l - Ax), (5) 

where A(x) and B(x) are analytic in the neighbourhood of x c . Note that the terms 
associated with the exponents 1/2 and 3/2 become part of the analytic correction to the 
(1 — 4x)~ 1//2 term. Near the singularity on the negative x-axis, x = x_ = — 1/4 we expect 
the singular behaviour 

T(x) ~C(x)(l + 4x) 13/2 , (6) 

where again C(x) is analytic near x_. We expect similar behaviour near the pair of 
singularities x = ±i/2, and finally at the roots of 1 + x + 7x 2 we expect the behaviour 
T(x) ~ D(x)(l + x + 7x 2 ) 2 log(l + x + lx 2 ). 

Next we turn our attention to the asymptotic behaviour of the coefficients of T(x). To 
standardise our analysis, we assume that the critical point is at 1. The growth constant of 
staircase and imperfect staircase polygons is 4, so we normalise the series by considering 
a new series with coefficients r n , defined by r n = t n+2 /4: n . Thus the generating function 
we study is 7Z(y) = J2 n>0 r n y n = 4 + 3y + 2.625y 2 + • • •. From equations © and © it 
follows that the asymptotic form of the coefficients is 

m) = r . * E + (_!,. (_«.)) + (A-»). (7) 

V i>0 ^ J 

The last term includes the effect of other singularities, further from the origin than the 
dominant singularities. These will decay exponentially since A > 1 in the scaled variable 
y = x/4. 

Using the recurrence relations for t n (derived from the ODE) it is easy and fast to 
generate many more terms r n . We generated the first 100000 terms and saved them as 
floats with 500 digit accuracy (this calculation took less than 15 minutes). With such a 
long series it is possible to obtain accurate numerical estimates of the first 20 amplitudes 
di, bi, Ci for i < 19 with precision of more than 100 digits for the dominant amplitudes, 
shrinking to 10-20 digits for the the case when i = 18, or 19. In making these estimates 
we have ignored the exponentially decaying term, which is the last term in eq.([Hj) . In 
this way we confirmed an earlier conjecture j3] that a = where we have taken 

into account the different normalisation, as discussed in the introduction. We also find 
that b = 3.173275384589898481765. . . and c = , though we have not been able to 
identify bo- However, we have successfully identified further sub-dominant amplitudes, 

nnrl finH n - ~ 89 n - 1019 nnr | n - -10484935 j _ 225 _ -16575 

and tind a x - ^37^ a 2 - ^75-^, and a 3 - 248832 ^ 3/ 2 > and c x - c 2 - j^^, 

and C3 = ^8^3 f 2 • ^ seems possible that the amplitudes 7r 3//2 v / 3aj and 7t 3//2 q are rational. 

Estimates for the amplitudes were obtained by fitting r n to the form given above using 
an increasing number of amplitudes. 'Experimentally' we find we need about the same 
total number of terms at x c and — x c = x_. 

So in the fits we used the terms with amplitudes a,, and bi, i = 0, . . . , K and q, 
i = 0, . . . , 2K. Going only to i — K with the c, amplitudes results in much poorer 
convergence and going beyond 2K leads to no improvement. For a given K we thus have 
to estimate 4K + 3 unknown amplitudes. So we use the last 4K + 3 terms r n with n 
ranging from 100000 to 100000 — AK — 2 and solve the resulting 4K + 3 system of linear 
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equations. We find that the amplitudes are fairly stable up to around 2K/3. We observed 
this by doing the calculation with K = 30 and K = 40 and then looking at the difference 
in the amplitude estimates. For ao and 60 the difference is less than 10~ 131 , while for cq 
the difference is less than 10~ 123 . Each time we increase the amplitude index by 1 we 
loose around 10 6 in accuracy. With i = 20 the differences are respectively around 10 -16 
and 10" 8 . 

The excellent convergence is solid evidence (though naturally not a proof) that the 
assumptions leading to equation (JJJ) are correct. Further evidence was obtained as follows. 
We can add extra terms to the asymptotic form and check what happens to the amplitudes 
of the new terms. If the amplitudes are very small it is highly likely that the terms are not 
truly present (if the calculation could be done exactly these amplitudes would be zero). 
One possibility is that our assumption about integer exponents leading only to analytic 
terms is incorrect. To test this we fitted to the form 

^E(^^ + M)"(^)) + o,A-»), 

(as above, in making these estimates we have ignored the exponentially decaying term, 
which is the last term in the above equation.) With K = 30 we found that the amplitudes 
d\ and 61 of the terms \ogn/n and 1/n, respectively, were less than 10~ 60 , while the 
amplitudes 03 and 63 were less than 10~ 50 . We think we can safely say that all the 
additional terms we just added are not present. We found similar results if we added 
terms like log 2 n or additional logra terms at y = —1. That is, we found that those 
terms were not present. So this fitting procedure provides convincing evidence that the 
asymptotic form (J7J), and thus the assumption leading to this formula, is correct. 

4 Conclusion 

We have developed an improved algorithm for enumerating three choice polygons. The 
extended series, coupled with a search program that assumes the solution is a Fuchsian 
ODE, enabled us to discover the underlying ODE, which was of 10 th order. We did this 
without using more than 50 of the coefficients that we had generated. That is to say, 50 
of the known coefficients were unused, and so their value provided a check on the solution 
found. This leads us to believe that we have found the correct ODE, as it reproduces 
the known, unused coefficients, though we have not provided a proof. Further refinement 
allowed us to find an 8 th order ODE. 

A numerical technique we have developed specifically for such problems then allowed 
us to find accurate numerical estimates for the amplitudes of the first several terms in the 
asymptotic form for the coefficients. 

We have also initiated an investigation of the area generating function. We expect this 
to involve g-series, and thus far our investigations only lead us to believe that the area 
generating function A(q) is of the form 

A(q) = (G(q) + H(q)/y/T=jfi)/[J (l, 1, q)\ 

where Jo(x,y,q) is a g-generalisation of the Bessel function, and occurs, for example, in 
the solution of the problem of staircase polygons enumerated by area p. Here q = 77 is 
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the first zero of Jq(1, 1, q), and G and H are regular in the neighbourhood of q = r\. The 
coefficients thus behave asymptotically as 

a n = [q n ]A(q) ~ const.i]~ n n 3/2 . 

In a subsequent publication [7] we propose to investigate the area generating function more 
fully, and to say more about the properties of the ODE we have found for the perimeter 
generating function. In particular, we discuss some simple solutions of the ODE, and ask 
what these can tell us about the full solution. 
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Appendix A: Polynomials Q n {x) 

Qs(x) = -180 + 4005a; - 45340a; 2 + 352567a; 3 - 2100653a; 4 + 8247059a; 5 + 1869782x 
-198745492a; 7 + 222232422a; 8 + 7981490552a; 9 - 58454247760a; 10 
+223070561538X 11 - 653903984242a; 12 + 1691567153918a; 13 
-3628069390936a; 14 + 9508812403200a; 15 - 42130737708796a; 16 
+151950842991736x 17 - 347187650580720a; 18 + 558723092175488a; 19 
-722483977609792a; 20 + 551434913787008a; 21 + 82817126361856a; 22 
-426478334005248x 23 + 279157576126464a; 24 + 2780644737024a; 25 



Q 7 (x) = -3420 + 82530a; - 926615a; 2 + 6866662a; 3 - 37878392a; 4 + 131975108a; 5 
+198512462x 6 - 5322566116a; 7 + 16816064102a; 8 + 88956629348a; 9 
-872972184658a; 10 + 3395585125316a; 11 - 8662194926872a; 12 
+2179593948608x 13 + 130585482759744a; 14 - 698610495175368a; 15 
+2229946022661696x 16 - 6216128747042864a; 17 
+ 15724091332879132a; 18 - 38607908490402392a; 19 
+ 128963713249678592a; 20 - 464640056155209952a; 21 
+1296873363475699328x 22 - 2966555758830491904a; 23 
+5741739615453110784a; 24 - 7824348079140616704a; 25 
+8096625038421797888a; 26 - 6327622359115208704a; 27 
-663175049190105088a; 28 + 4390942020748738560a; 29 
-3449431865352388608a; 30 - 33011814317948928a; 31 
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Q 6 (x) = -9180 + 310275a; - 4493475a; 2 + 40204094a; 3 - 262917778a; 4 

+ 1302960911a; 5 - 3743237840a; 6 - 8573351756a; 7 + 140454430666a; i 
-409322626730a; 9 - 1570504457342a; 10 + 18303308342032a; 11 
-89658228463172a; 12 + 259420736216632a; 13 + 26862202296376a; 14 
-4190021721023184a; 15 + 21897720821926584a; 16 
-75837533674259508a; 17 + 212508813586272428a; 18 
-476010656497826944a; 19 + 1034090705056496672a; 20 
-3196181326637410304a; 21 + 10833216991064882848a; 22 
-30172750280212408832a; 23 + 70340668591569812736a; 24 
-132300506186507025408a; 25 + 177280513560453634560a; 26 
-184990945124657242112a; 27 + 135828858351882342400a; 28 
+ 12754320650381836288a; 29 - 85576383794502107136a; 30 
+61165017902554546176x 31 + 565560894352785408a; 32 



Q 5 (x) = -4500 + 244800a; - 4876845a; 2 + 55164150a; 3 - 438701640a; 4 

+2758453094a; 5 - 13804842198a; 6 + 45370091528a; 7 - 3608230380a; 8 
-892524490064a; 9 + 4421327158154a; 10 - 2297315126532a; 11 
-103201897035096a; 12 + 748998082407080a; 13 - 2329708885595260a; 
-457382726191024a; 15 + 35817660448173240a; 16 
-188156345496838984a; 17 + 677783573996257364a; 18 
-1904649390940935752a; 19 + 4199594693024922016a; 20 
-8814226144821806432a; 21 + 23568486792872894272a; 22 
-70089404940793421632a; 23 + 188311273940137111552a; 24 
-435002993494719438848a; 25 + 791152555777632593920a; 26 
-1045593345640931730432a; 27 + 1096015208846337957888a; 28 
-774016903940080771072a; 29 - 37178029375778357248a; 30 
+412071049964952354816a; 31 - 275345921075326746624a; 32 
-2464051649845395456a; 33 



Q 4 (x) = 31500 - 1114080a; + 17560755a; 2 - 178469565a; 3 + 1412918104a; 4 

-9431590849a; 5 + 52336335969a; 6 - 220707961458a; 7 + 52596571 1332a; ; 
+915935968370a; 9 - 13996439933349a; 10 + 35303141246088a; 11 
+231992664240696a; 12 - 2180352456480752a; 13 + 6859298731027888a; lz 
-1272161420555012a; 15 - 75338205421491734a; 16 
+406836568590013948a; 17 - 1513874477368697252a; 18 
+4439738234446975124a; 19 - 10514406278248398472a; 20 
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+22797584086521040520x 21 - 52624647215757093584a; 22 
+130673617185226821792x 23 - 324680301683722155712a; 24 
+724969297042825531136a; 25 - 1271869215082051692800a; 26 
+1661614177465373698560a; 27 - 1744486537247742479360a; 28 
+1209724637255295010816a; 29 - 32142663152460406784a; 30 
-498040622799430975488a; 31 + 321317702703841148928a; 32 
+2787318284392857600x 33 



Q 3 (x) = -156000 + 3778920a; - 46727325a; 2 + 457371630a; 3 - 3919246431a; 4 

+27446185200a; 5 - 152613919692a; 6 + 659637747242a; 7 - 1723470963068a; 6 
-1667066145852a; 9 + 27889854017778a; 10 + 15933308039400a; 11 
-972460279627326a; 12 + 4552136023731292a; 13 - 7976188460233924a; 14 
-4422880527966948a; 15 + 63325989574562728a; 16 
-287206984863975352a; 17 + 1115308575007981980a; 18 
-3508943115779966584a; 19 + 8987842561562515768a; 20 
-19184807012355087408a; 21 + 37821550927408731776a; 22 
-83609060792238083072a; 23 + 194683017390969665280a; 24 
-366258975512082319872a; 25 + 498254429378056694784a; 26 
-558421919820222289920a; 27 + 441211762632912959488a; 28 
-80404063142199537664a; 29 - 110342796490113417216a; 30 
+85904442856027127808x 31 + 720965567415582720a; 32 



Q 2 (x) = 102000 - 1245240a; + 445275a; 2 + 77507430a; 3 - 505005638a; 4 + 674357270a; 5 
+ 7410398802a; 6 - 50751541108a; 7 + 109730141494a; 8 + 263567061768a; 9 
-2398666258514a; 10 + 4447124418524a; 11 + 33544348232760a; 12 
-341405641395740a; 13 + 1843130781900080a; 14 - 7441271357292384a; 15 
+23827305830694324a; 16 - 59142500096057112a; 17 
+ 113845825936073424a; 18 - 169659492965796928a; 19 
+ 190085091 157739584a; 20 - 160391840217609984a; 21 
+95477320250924800a; 22 - 21461546279272960a; 23 
-73590898428536832a; 24 + 43442402559821824a; 25 
+ 129164030193680384a; 26 - 136460131311329280a; 27 
+5453275269051 1872a; 28 + 389290263183360a; 29 



Qo(x) = Qi{x) = -Q 2 {x) 
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